load c:\manoj\projects\plasma\HUA_PIU HUA_PIU_fday HUAPIU magnetic_frac;
load c:\manoj\projects\plasma\Julia_PLP.mat;
bin_radius = 30/(60*24);%<------- NORMAL
%bin_radius = 60/(60*24);
kp_limit = 3;
lt_step =  1;
%for delay_index = 1:13,
delay_index = 2;
delay_sel = [10,13,14,15,16,17,18,19,20,21,22,23,24];
Kp_delay = delay_sel(delay_index);%Valid columns [10,13-24];
kp_delay_text = ['00';'00';'01';'02';'03';'04';'05';'06';'07';'08';'09';'22';'MN'];
lag_step = 30;
min_data = 10; %Minimum number of data pairs to be considered for correlation
%LAG = [0:lag_step:360]/(60*24);
LAG = [30:lag_step:240]/(60*24);%<------- NORMAL
LT = [10:lt_step:17]; %<------- NORMAL
%LT = [9:lt_step:17];
plot_option = 2;%[1 scatter, 2 image 3 curves];
ED_matrix = ones([length(LT)-1,length(LAG)])*NaN;
LV_matrix = ones([length(LT)-1,length(LAG)])*NaN;
RT_matrix = ones([length(LT)-1,length(LAG)])*NaN;
PT_matrix = ones([length(LT)-1,length(LAG)])*NaN;
ED_error  = ones([length(LT)-1,length(LAG)])*NaN;
LV_error  = ones([length(LT)-1,length(LAG)])*NaN;
RT_error  = ones([length(LT)-1,length(LAG)])*NaN;
PT_error  = ones([length(LT)-1,length(LAG)])*NaN;

N_good = 1;
for jj = 1:length(LAG),
for kk = 1:length(LT)-1,
    
     L = ~isnan(PLP(:,27)) & ~isnan(PLP(:,18))  &  D_KP(:,Kp_delay) <=kp_limit & PLP(:,7) >= LT(kk) & PLP(:,7) <= LT(kk+1)...
         & PLP(:,12) == 1;% & PLP(:,20) >= 1.5;%& PLP(:,29) >= 1.5;% & PLP(:,8) > -85 & PLP(:,8) < -70;
  %  L = ~isnan(PLP(:,18))  & PLP(:,10) <=kp_limit & PLP(:,7) >= LT(kk) & PLP(:,7) <= LT(kk+1) & PLP(:,12) == 1;% & PLP(:,8) > -85 & PLP(:,8) < -70;
         IndData = PLP(L,:);
         counter = 0;
         
         for i = 1:size(IndData,1),
         fday_champ = tmjd(IndData(i,1),IndData(i,2),IndData(i,3),IndData(i,4),IndData(i,5),IndData(i,6),1);   
         iday = find(HUA_PIU_fday==floor(fday_champ));
         if ~isempty(iday),
            fday_champ_frac =  fday_champ - floor(fday_champ) - LAG(jj);
            LLL = isnan(HUAPIU(iday,:));kkk = find(LLL == 0);
            if length(kkk) > 3,
            if fday_champ_frac > magnetic_frac(kkk(1)) & fday_champ_frac < magnetic_frac(kkk(end)),
                 counter=counter+1;
                [trash left_bin] =  min(abs(magnetic_frac(kkk(1):kkk(end)) - (fday_champ_frac-bin_radius)));
                [trash right_bin] = min(abs(magnetic_frac(kkk(1):kkk(end)) - (fday_champ_frac+bin_radius)));
                left_bin  = kkk(1) + left_bin - 1;
                right_bin = kkk(1) + right_bin -1;
               HUA_data(counter) = nanmean(HUAPIU(iday,left_bin:right_bin));%<------- NORMAL
                %HUA_data(counter) = nanmean(HUAPIU(iday,800:1300)); %
                Plp_ed(counter) = (IndData(i,18) + IndData(i,27))/2;
                %Plp_ed(counter) = (IndData(i,18)); %WATCH OUT ONLY ONE PEAK !
                Plp_lv(counter) = (IndData(i,17) + IndData(i,26))/2;
                %Plp_ed(counter) = (IndData(i,18) + IndData(i,27))/2;
                Plp_ratio(counter) = (IndData(i,20) + IndData(i,29))/2;
                Plp_pt(counter) = Plp_ed(counter) -IndData(i,36);
                Plp_altitude(counter) = IndData(i,9);
                plot(magnetic_frac,HUAPIU(iday,:),magnetic_frac(left_bin:right_bin),HUAPIU(iday,left_bin:right_bin));
                title(sprintf('fday %8.3f lt %5.2f long %6.2f CTR %4.2f',...
                    fday_champ,IndData(i,7),IndData(i,8),Plp_ratio(counter)));
                pause;
               
            end;
            end;
         end;
         end;
            
         lt_count(kk,jj) = counter;
        % fprintf('LT=%5d LAG = %10.5f COUNT = %3d',LT(kk),LAG(jj)*(60*24),counter); 
         if counter >= 1,
             [cof_ed,err_ed] = corrcoef(HUA_data,Plp_ed); 
             [cof_lv,err_lv] = corrcoef(HUA_data,Plp_lv); 
             [cof_rt,err_rt] = corrcoef(HUA_data,Plp_ratio); 
             [cof_pt,err_pt] = corrcoef(HUA_data,Plp_pt); 
             ED_matrix(kk,jj) = cof_ed(2,1);
             LV_matrix(kk,jj) = cof_lv(2,1);
             ED_error(kk,jj) = err_ed(2,1);
             RT_matrix(kk,jj) = cof_rt(2,1);
             RT_error(kk,jj) = err_rt(2,1);
             PT_matrix(kk,jj) = cof_pt(2,1);
             PT_error(kk,jj) = err_pt(2,1);
             LV_error(kk,jj) = err_lv(2,1);
             
             if plot_option == 1,
             if cof_rt(2,1) < 100,
             subplot(8,8,(jj-1)*(length(LT)-1)+kk);
             %set(gca,'FontSize',16);
             hold on;
             plot(HUA_data,Plp_ratio,'r.','MarkerSize',15);
             grid on;
             axis([-10 100 1 3]);
             text(30, 1.2,sprintf('%5d',LT(kk)));
             text(30, 1.8,(max(Plp_ratio)-min(Plp_ratio))*0.3+min(Plp_ratio),sprintf('%5d',LAG(jj)*(60*24)));
             text(30, 2.4,(max(Plp_ratio)-min(Plp_ratio))*0.7+min(Plp_ratio),sprintf('%4.2f',cof_rt(2,1)),'color','r');
             %xlabel(sprintf('%5d-%5d',LT(kk),LAG(jj)*(60*24)));
             N_good = N_good+1;
             end;
             end;
        end;
%         fprintf('\n');
         clear HUA_data Plp_ed Plp_lv Plp_ratio Plp_pt;
      
end;
end;

for i = 1:length(LT)-1,
SS(i,:) = sprintf('%2d-%2d',LT(i),LT(i+1));
end;

if plot_option == 3,
color_index = ['k*-'; 'b*-' ;'r*-'; 'g*-'; 'c*-' ;'y.-' ;'m.-'; 'r>-' ;'g.-'; 'c.-'];

% color_map = colormap('jet');
% color_index = color_map([1:3:64],:);
figure(1);
subplot(331);
for i = 1:length(LT)-1,
    hold on;
errorbar(LAG*(24*60),PT_matrix(i,:)',PT_error(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
xlabel('LAG in minutes');
ylabel('CC');
hold on;
axis([LAG(1)*24*60 LAG(end)*24*60 0 1]);
end;
title('C-T - Julia correlation');
text(-42,1.0,['Kp-' sprintf('%d',kp_limit)]);
text(-42,0.8,['Delay-' kp_delay_text(delay_index,:)]);
subplot(334);
for i = 1:length(LT)-1,
plot(LAG*(24*60),lt_count(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);

hold on;
end;
xlabel('LAG in minutes');
ylabel('Number of data points');

subplot(337);
for i = 1:length(LT)-1,
plot(LAG*(24*60),log10(PT_error(i,:))','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Log |Error|');


subplot(332);
for i = 1:length(LT)-1,
    hold on;
errorbar(LAG*(24*60),RT_matrix(i,:)',RT_error(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
xlabel('LAG in minutes');
ylabel('CC');
hold on;
axis([LAG(1)*24*60 LAG(end)*24*60 0 1]);
end;
title('CTR - Julia correlation');

subplot(335);
for i = 1:length(LT)-1,
plot(LAG*(24*60),lt_count(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);

hold on;
end;
xlabel('LAG in minutes');
ylabel('Number of data points');

subplot(338);
for i = 1:length(LT)-1,
plot(LAG*(24*60),log10(RT_error(i,:))','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Log |Error|');


subplot(333);
for i = 1:length(LT)-1,
    hold on;
errorbar(LAG*(24*60),LV_matrix(i,:)',LV_error(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
xlabel('LAG in minutes');
ylabel('CC');
hold on;
axis([LAG(1)*24*60 LAG(end)*24*60 0 1]);
end;
title('L-Value - Julia correlation');

subplot(336);
for i = 1:length(LT)-1,
plot(LAG*(24*60),lt_count(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);

hold on;
end;
xlabel('LAG in minutes');
ylabel('Number of data points');

subplot(339);
for i = 1:length(LT)-1,
plot(LAG*(24*60),log10(LV_error(i,:))','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Log |Error|');

end;